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1 .        INTRODUCTION 

We  have  investigated  several  time-frequency  decomposition  techniques  and  their 
applications  in  the  detection  and  analysis  of  digitally  modulated  signals  [10],  [14],  [13].  In  this 
report,  we  will  specifically  discuss  wavelet  decompositions  and  present  a  few  preliminary  results. 

The  one  particular  digital  modulation  that  we  consider  in  here  is  binary  phase-shift-key 
(BPSK)  signal,  which  is  a  sinusoid  modulated  by  a  baseband  rectangular  data  waveform.  Over  a 
duration  of  one  binary  data  bit,  the  phase  of  a  BPSK  wave  is  either  0°  or  180°  depending  on  the 
value  of  the  data,  which  is  either  +1  or  -1.  For  example,  a  wave  with  data  rate  9600  bps  can 
change  phase  every  1/9600  seconds.  In  tactical  communications  using  direct-sequence  spread- 
spectrum  (DS/SS)  [1 1],  the  corresponding  spread  BPSK  wave  can  change  phase  N  times  faster, 
where  N  is  the  spectrum  spreading  factor.    Since  BPSK  and  its  extensions  (such  as  QPSK  and 
OQPSK)  are  the  most  predominant  modulation  formats  used  in  practice,  there  are  interests  in 
developing  efficient  methods  for  detecting  the  phase  shifts  under  various  channel  and  noise 
scenarios.  Previous  research  efforts  have  used  cyclostationarity  [6]  and  ambiguity  function  with 
specific  kernel  as  the  smoothing  function  [23],  [13]. 

In  Chapter  2,  we  will  discuss  the  so-called  continuous  wavelet  transform  (C WT).  In 
Chapter  3,  we  will  discuss  the  discrete  wavelet  transform  (DWT)  and  the  computations  of 
orthonormal  compactly  supported  scaling  and  wavelet  functions  on  the  dyadic  numbers.  In 
Chapter  4,  we  use  the  wavelet  time-frequency  decomposition  techniques  to  detect  ±180°  phase 
shifts  in  noise.  A  few  concluding  remarks  are  given  in  Chapter  5. 


2 .        CONTINUOUS  WAVELET  TRANSFORM 

The  continuous  wavelet  transform  (CWT)  of  a  L2(E)  signal  x(t)  is  defined  as  [8],  [2],  [3]: 


».  d-it^  *<=*)* 


(2.1) 


where  y(t)  is  the  wavelet  function, "  *  "  denotes  complex  conjugation,  and  a  >  0.  The  scaling 
factor  "a"  in  (2.1)  is  taken  to  be  power  of  2  in  this  report,  i.e.,  a  =  2_m,  for  m  =  0,  1,  2,  3  ,  .  . 
The  inverse  transform  is  given  by: 
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where  the  constant 

=   f  i  *vm  i2 
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must  be  finite.  This  implies  that  the  wavelet  function  must  satisfy 

¥(f)  =  0,        i.e.,    J  y(t)  dt  =  0  (2.4) 

Thus,  the  wavelet  function  is  bandpass  type  with  Fourier  transform  ^(f)  =  0  at  f  =  0  ; 
equivalently,  the  area  under  the  wavelet  function  is  exactly  0.  Note  that  we  can  rewrite  (2.1)  as 

X(t,  a)  =  Va  J   X(f)  4""(af)  exp(j27tft)  df  (2.5) 

in  terms  of  the  Fourier  transforms.  The  CWT  can  be  computed  via  (2.1)  or  (2.5). 

The  CWT  X(t,  a)  represents  a  decomposition  of  x(t)  into  a  time-scale  plane.  Usually, 
X(t,a)  is  graphed  against  time  t  and  the  log  scale  log2(a).  The  CWT  can  be  regarded  as  a  constant- 
Q  bandpass  filter  since  (2. 1 )  is  a  convolution  of  two  time  functions.  The  linear  bandpass  filter 
impulse  response  at  level  m  is 

\|/m(t)  =  2m/2  y*(2mt)  (2.6) 

If  we  define  "fm"  as  the  r.m.s.  center  frequency  of  the  passband  by 

oo 

f    (f-fm)2l^m(f)l2df  =0 

Jo  (2.7) 


then  we  have 

fm  =  2mf0      and     am  =  2ma0  (2.8) 

where  fo  is  the  r.m.s.  center  frequency  and  2an  is  the  r.m.s.  bandpass  bandwidth  of  the  basic 
wavelet  function  \jf(t).  The  ratio  of  fm  to  am  is  always  equal  to 

am      ao  (2.9) 

which  is  constant  for  all  decomposition  level  m.    As  the  level  of  decomposition  (i.e.,  m)  increases, 
the  filtering  function  \|/m(t)  of  (2.6)  becomes  narrower  in  time.  The  time-resolution  of  the  wavelet 
decomposition  therefore  increases  with  m.  On  the  other  hand,  the  r.m.s.  bandpass  bandwidth 
(i.e.,  2am  )  of  \|/m(t)  increases  exponentially  with  m,  and  therefore  the  frequency-resolution  of  the 
wavelet  decomposition  is  better  for  smaller  m.    Unlike  the  classical  time-frequency  transforms 
such  as  the  short-time  Fourier  transform  (STFT)  or  the  Gabor  transform  (GT)  [1],  [9],  which  do 
not  scale  the  filtering  function  and  and  have  constant  ume-resolution  at  all  frequencies,  the  wavelet 
transform  has  higher  time-resolution  at  higher  frequency.  One  can  observe  that  for  CWT,  its 
time-resolution  is  increasingly  better  at  higher  frequencies,  while  its  frequency-resolution  is 
increasingly  better  at  lower  frequencies. 

We  will  use  the  complex  modulated  Gaussian  wavelet  function: 

\|f(t)  =  exp(j27ifVj/t)  exp(-t2/2a2) 

=  cos  2?cfvt  exp(-t2/2a2)  +  j  sin  2rcfvt  exp(-t2/2a2)  (2.10) 

Of  course  the  Fourier  transform  of  a  Gaussian  function  is  a  Gaussian  function,  and  the  complex 
exponential  time  function  exp(j27rfyt)  will  cause  a  shift  by  fy  in  frequency  in  the  Fourier 
transform.  It  is  then  possible  to  choose  fy  big  enough  so  that  for  complex  modulated  Gaussian 
wavelet  we  have 

4/(f)  =  0      for   f  <0  (2.11) 

Condition  (2.1 1)  is  especially  desirable  when  the  wavelet  is  used  to  analyze  sinusoids.  For 
example,  let  x(t)  =  cos  27tfct,  for  t  in  a  very  large,  yet  finite,  window.  Using  (2.5),  we  have 

X(t,  a)  =  Y  exp(j27tfct)  ¥*(afc)  +  ^  exp(-j27tfct)  ¥*(-afc)  (2. 1 2) 


If  4^0  satisfy  (2. 1 1),  then  the  second  term  in  (2. 12),  which  is  a  contribution  by  the  impulse 
5(f+fc),  drops  out.  It  follows  that  the  phase  of  the  resulting  CWT  X(t,  a)  will  nicely  reveal  the 
center  frequency  of  the  carrier  wave  x(t)  =  cos  27tfct.    For  other  stationary  or  non-stationary 
signals,  magnitude  and  phase  of  X(t,  a),  with  \y(t)  being  the  complex  Gaussian  wavelet,  will  reveal 
the  frequency  as  well  as  the  time  behaviors  of  the  signal. 

We  also  remark  that  we  can  multiply  a  suitable  baseband  function  b(t)  with  the  complex 
sinusoid  exp(j27tfyt)  to  obtain  a  wavelet  function: 

y(t)  =  exp(j27ifvt)  b(t)  (2. 1 3) 

so  that  its  Fourier  transform  is  zero  on  the  negative  frequency  axis.  The  above  construction  yields 
a  large  class  of  potentially  useful  wavelet  functions.  The  preference  to  use  the  Gaussian  envelope 
exp(-t2/2a2)  is  due  to  the  fact  that  the  Gaussian  function  achieves  the  lower  bound  in  the 
uncertainty  principle: 

AtAf>—  (2.14) 

47t 

with  equality,  where  At  is  the  r.m.s.  time-spread  and  Af  is  the  r.m.s.  frequency-spread  of  the 
signal  [20].    Note  that  two  pulses  in  time  can  be  discriminated  only  if  they  are  more  than  At 
seconds  apart  in  time,  and  two  sinusoids  can  be  discriminated  only  if  they  are  more  than  Af  apart 
in  frequency.  Equation  (2.14)  states  that  one  can  only  trade  time-resolution  for  frequency- 
resolution,  or  vice  versa. 

Using  CWT,  the  frequency-resolution  Af  is  in  the  order  of  am  =  2m  Gq.  Viewing  the 
CWT  X(t,  a)  in  a  time-scale  plane  allow  us  to  look  at  the  signal  at  different  scale  (i.e.,  different 
frequencies  in  the  passband  of  2am  Hz  wide  around  the  center  frequency  fm).  The  CWT  indeed 
presents  a  multiresolutional  analysis  of  the  signal,  as  if  the  signal  is  being  processed  by  a  bank  of 
filters  with  logarithmically  spaced  centered  frequency  fm  but  constant  bandwidth  to  frequency  ratio 
of  20m/fm  =  2ao/fo-  These  relationships  are  best  visualized  by  considering  the  complex  modulated 
sine  wavelet  function: 

\|/(t)  =  exp(j2rcfM/t)  -=■  sine  ^  (2. 1 5) 


which  has  a  rectangular  passband  with  bandwidth  equal  to  2gq  =  W.  The  center  frequency  of  \j/(t) 
is  exactly  equal  to  fy.  We  also  propose  to  use  this  complex  modulated  sine  wavelet  to  analyze  our 
signals. 

Computation  of  X(t,  a)  could  be  lengthy  and  expensive.  In  practice,  however,  we  only 

have  a  finite-length  sampled  sequence:  x(nT),  n  =  0,  1,  2, ,  Ns-1,  where  Ns  is  the 

sequence  length  and  the  sampling  rate  is  at  least  equal  to  the  Nyquist  rate,  which  is  twice  the 
highest  frequency  component  in  the  signal  x(t).  We  will  only  compute  X(t,  a)  at  the  sampling 
times  t  =  nT  and  at  levels  m  =  0,  1,  2, .  . . . ,  [log2  Ns].  The  integral  of  (2. 1)  is  now  reduced  to  a 
finite  summation  and  the  computation  burden  is  greatly  relieved.  In  this  discrete-time 
approximation  of  the  CWT,  we  essentially  perform  a  multiresolutional  analysis  of  an  auxiliary 
signal  xa(t),  which  is  a  piecewise  linear  function  constructed  by  joining  two  consecutive  sampled 
points  x(nT)  and  x((n+l)T)  with  a  straight  line.  Towards  this  end,  the  CWT  of  (2.1)  produces  a 
total  of  Nmax  =  [log2  Ns]  +  1  snap-shots  of  the  auxiliary  signal  via  bandpass  filterings. 

The  complex  modulated  Gaussian  wavelet  of  (2. 10)  is  a  nice  analytical  function.  It  is 
continuous  and  differentiable  everywhere.  Perhaps  one  drawback  is  that  it  has  an  infinitely  span 
and  therefore  it  is  not  a  causal  filter  impulse  response.  In  practice,  however,  this  mishap  can  be 
resolved  by  using  a  wavelet  with  sufficiently  fast  decay  over  the  processing  window  so  that  it  can 
be  regarded  as  of  finite  span.   We  select  the  Gaussian  wavelet  because  of  two  additional  reasons: 
(a)  tunability;  and  (b)  bandwidth  adjustability.  By  varying  the  modulating  frequency  fv,  we  can 
place  the  passband  anywhere  on  the  frequency  axis.  Second,  by  altering  a,  the  bandwidth  of  the 
wavelet  function  \j/(t)  can  be  selected  at  will.  Each  (fy,  a)  selection  will  result  in  a 
multiresolutional  analysis  that  consists  of  Nmax  snap-shots  of  the  signal.    A  systematic  selection 
of  (fy,  a)  will  allow  a  thorough  analysis  of  the  signal  at  all  possible  time-resolution  and 
frequency-resolution  levels. 

We  show  an  example  of  a  complex  modulated  Gaussian  wavelet  function  \j/(t)  with  fy  =  8 
Hz  and  a  =  0. 12  in  Figure  2. 1 .  This  value  is  selected  so  that  all  the  energy  in  the  wavelet  is 
essentially  contained  in  the  compact  time  interval  [-0.5, 0.5]  for  presentation  purpose.    In  the  left 
column  of  Figure  2. 1 ,  the  real  part,  imaginary  part,  magnitude,  and  phase  of  the  wavelet  are 
depicted.  In  the  right  column,  the  energy  spectral  densities  (which  is  magnitude  square  of  the 
Fourier  transform)  of  the  real  part,  imaginary  part,  and  the  wavelet  function  itself  are  shown.  Note 
that  this  complex  wavelet  has  no  spectral  component  in  the  negative  frequency  axis,  a  desirable 
property  as  discussed  previously.  The  spectrum  is  centered  around  f  =  fy  =  8  and  the  r.m.s. 
passband  bandwidth  of  \|/(t)  is  2c  -  0.24.    In  processing  a  given  signal  x(t),  the  signal  is  first 
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normalized  to  be  contained  in  the  time  interval  [0,  1];  it  is  then  decomposed  by  the  analyzing 
functions  ym(t),  with  \|/o(t)  =  \|/(t)  contained  also  in  the  interval  [0,  1]. 
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Fig.  2.1  Complex  modulated  Gaussian  wavelet  function  \j/(t):  time-domain  functions  and 

energy  spectral  densities,  fy  =  8,  a  =  0.12. 


10 
3.       DISCRETE  WAVELET  TRANSFORM 

In  this  section  we  will  discuss  an  orthonormal  decomposition  of  a  signal  x(t)  using  a 
scaling/wavelet  function  pair,  and  describe  the  so-called  discrete  wavelet  transform  (DWT) 
algorithm  using  orthonormal  wavelet  functions  [4],  [16],  [17],  [18].  We  will  also  discuss  the 
multiple-phase  DWT,  and  generations  of  wavelet  functions  that  have  compact  (i.e.,  closed  and 
bounded)  support. 

3.1      Orthonormal  Wavelet  Decomposition 

Consider  a  nested  sequence  of  closed  subspaces  (Vm;  me  2}  in  L2(E)  such  that 

{0}  =  -  •cV.2cV.,cV0cV,cV2c  •  =L2(R)  (3.1) 

where  O  Vm  =  {0}and    0    Vm  =  L  (R).  One  can  find  an  unique  scaling 

me  Z  me  Z 

(or  averaging)  function  <}>(t)  such  that  (<}>mtn(t);  ne  Z} 
is  an  orthonormal  basis  of  Vm  where 


WO  =  a^  *<  «mt  ~  nT)  (32) 

by  normalizing,  dilating  or  contracting,  and  time-shifting  <{>(t).  The  value  of  T  (>0)  is  arbitrary.  As 
is  common  in  dyadic  multiresolution  analysis,  we  take  a  (*0,1)  to  be  2  in  this  report. 

The  orthogonal  projection  of  a  finite-energy  deterministic  function,  x(t),  onto  the  subspace 
Vm  is  an  approximation  of  that  function  at  resolution  level  m.  We  can  write: 


neZ 

where 

bm,n  =  <  x(t),  +mJfi)  >  =  f    x(t)  Cn  (Odt    =  (  x  *  tm,0  )(n2"mT) 

J-~  (3.3) 

is  a  coefficient  equal  to  the  orthogonal  projection  of  x(t)  onto  the  basis  function  (j)m,n(t)-  Here,  x(t) 
=  x(-t).  Since  {Vm}  is  a  sequence  of  densed  and  nested  subspaces,  we  have  xm(t)  — >  x(t)  in 
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L2(R)  as  m  — »  <».  By  dilating  (when  m<0)  or  contracting  (m>0)  the  scaling  function, 
approximation  for  x(t)  at  any  resolution  level  can  be  obtained  by  increasing  the  number  of 
coefficients  exponentially  with  m.   Although  the  set  of  functions  {<j>m,n(t);  (m» n)  e  22}  spans 
the  entire  space  L2(IR),  it  is  not  a  minimal  set  because  the  subspaces  { Vm}  are  nested;  in 
particular,  (J^nCO  and  <J>i,k(t)  are  correlated  for  m  *  X.  For  many  applications,  for  example  in 
image  processing  and  compact  coding,  it  is  desirable  to  construct  an  orthonormal  multiresolution 
basis  for  L2(R).  Towards  this  end,  let  Vm  =  Vm.i  ©  Wm.i,  where  Wm.i  is  the  orthogonal 
complement  of  Vm_i  in  Vm.  Consider  Vq  with  {<J>(t-nT);  n  e  2 }  as  its  orthonormal  basis.  One 
can  construct  a  wavelet  function  \y(t)  from  <J>(t)  so  that  {\j/(t-nT);  n  e  Z }  is  an  orthonormal  basis 
for  Wo;  furthermore,  { \J/m,n(0;  n  e  2 }  =  { 2m/2\|/(2mt  -  nT);  n  e  Z }  is  an  orthonormal  basis  for 
Wm,  and  the  entire  set  {\ym,n(0;  (m,n)  e  Z2}  is  an  orthonormal  multiresolution  basis  for  L2(IR). 
The  relationship  between  the  subspaces  Vs  and  Ws  is  depicted  in  Fig.  3.1. 
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Figure  3.1        Subspaces  Vs  and  Ws. 


The  orthogonal  projection  of  x(t)  onto  the  subspace  Wm  is  the  detail  function  at  resolution 


level  m: 


dm(t)  =  ^dmtn\|/m,n(t), 

neZ 

=  xm+1(t)  -xm(t), 


dm,n  =  <x(t),\|/m,n(t)> 


(3.4) 
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This  detail  function  constitutes  the  difference  between  the  approximation  function  of  x(t)  at 
resolution  level  m+1  and  resolution  level  m.  Thus,  the  coefficients  {bm,n;  n  e  2 }  capture  the 
information  about  x(t)  at  resolution  level  m;  while  the  coefficients  (dm)n;  n  e  Z }  capture  the  detail 

information  going  from  resolution  level  m  to  resolution  level  m+1.  Note  that  L2(IR)  =  mm^  Wm 

and  the  Ws  are  mutually  orthogonal  subspaces.  The  sequences  {bm,n;  n  e  Z }  and  {dm)I1;  n  e  2 } 
can  be  obtained  directly  through  linear  filtering  of  x(t)  and  then  sample  at  time  t  =  n2_mT  (Fig. 
3.2).  Note  that  ^[(j)     (-t)]  yields  a  lowpass  transfer  function,  while  &[\\f    Q(-t)J  is  bandpass. 

* 
The  wavelet  decomposition  filters  \\f     (-t)  can  be  interpreted  as  a  parallel  bank  of  bandpass  filters 

(with  overlapping  passbands  for  adjacent  resolution  levels)  that  have  constant  size  passband  on  the 
logarithm  scale. 

t  =  n2"MT 


x(t) 


<W  H) 

#" 

VHO  <r0 

t  =  n2"mT 


"►    bM,n 


-►    dHn 


•►    4 


m,n 


Figure  3.2        Orthonormal  wavelet  decomposition  of  the  signal  x(t). 
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3.2      Discrete  Wavelet  Transform 

There  exist  recursive  (pyramidal)  algorithms  that  efficiently  compute  the  lower  resolution 
level  coefficients  {bi>n;  n  e  Z }  and  {djtiTi;  n  e  2 }  from  the  higher  resolution  level  coefficients 
{bm,n;  n  e  Z },  i  <  m;  conversely,  the  bs  at  the  higher  resolution  level  m  can  be  computed  from 
the  bs  at  a  base  (lowest)  resolution  level  M  and  all  the  ds  at  lower  resolution  levels  i,  M  <  i  <  m 
These  pyramidal  algorithms  are  the  so-called  discrete  wavelet  transform  (DWT)  and  they  are 
described  by  the  block  diagrams  in  Fig.  3.3. 


(a)  resolution 


m-l,n 


2x-downsample  decimators 
(keep  every  other  samples) 


(b)  reconstruction 


2x-upsample  interpolators 
(insert  "zero"  between  each  sample) 


Figure  3.3        Discrete  wavelet  transform. 


In  Fig.  3.3,  h(k)  is  the  impulse  response  of  a  discrete-time  conjugate  filter  that  satisfies 

IH(coT)  |2  +  I  H(coT+7t)  I2  =  1  (3.5) 


and  g(k)  is  the  impulse  response  of  the  mirror  filter  of  H  with 

G(co)  =  exp(-jw)  H*(o>+7r.)  (3.6a) 

and       H(G)T)G*((oT)  +  H(coT  +  7t)G*((oT  +  7t)  =  0  (3.6b) 

One  common  choice  of  g(n)  that  satisfies  (3.6)  is  g(n)  =  (-l)1-n  h(l-n).  The  scaling  function  and 
the  wavelet  function  are  orthogonal  to  each  other  and  they  are  related  by 

<D(2co)  =  H(co)  <I>(co) 

*P(2co)  =  G(CO)  <D(03)  (3.7) 

It  follows  that  O(co)  =  n~=]  H(2-m(o);  depending  on  the  choice  of  H(co),  it  is  possible  to  obtain 

scaling  function  <J>(t)  that  has  good  localization  properties  in  both  the  frequency  and  the  spatial 
domain.  The  properties  of  H  and  G  in  Mallat's  pyramidal  algorithms  can  be  derived  without 
reference  to  the  multiresolution  analysis  [2]. 

Knowing  the  approximation  function  xm(0  and  all  the  detail  functions  di(t)  at  resolution 
level  X  as  M,  M+l,  M+2,  .  .  .  ,  m-2,  m-1,  is  necessary  and  sufficient  to  construct  the 
approximation  function  xm(t).  The  approximation  of  x(t)  at  resolution  level  m  can  be  written  as 

Xm(t)  =  X  bm.n  <L.n<t)  =   X  bM,n  Wl>  +       I        I  di,n  Vi.„(t) 

neZ  neZ  M<i<mneZ  (3.8) 

for  any  M  <  m.  Without  loss  of  generality  (WLOG),  one  may  set  M  =  0  and  (J)(t)  =  4>o,o(t)  (tne 
lowest  resolution  function  physically  realizable),  and  consider  resolution  level  0  <  m  <  mmax-  In 
this  report,  we  leave  M  as  an  arbitrary  integer.  An  example  pair  for  ty  and  \\f  is  the  indicator 
function  on  [0,  T)  and  the  Haar  wavelet  on  [0,  T): 


<j)(t)  = 


f  i/Vt     ,  o  <  t  <  t 

V(t)  = 
0  ,  otherwise 


1/VT       ,  0  <  t  <  T/2 

(3.9) 
-WT      ,  T/2  <  t  <  T 
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In  general,    J  4>(t)  dt  =  <I>(0)  =  VT,  <J>(±°°)  =  0,  and    J  \j/(t)  dt  =  ¥(0)  =  0;  the  scaling 

-oo  -oo 

function  <|>(t)  can  be  interpreted  as  a  lowpass  function  while  \y(t)  is  bandpass.  The  Haar  wavelet 
forms  the  only  orthonormal  basis  of  compactly  supported  wavelets  for  which  the  associated 
scaling  function  <J>  has  a  symmetry  axis.  Also  note  that  the  Haar  wavelet  is  not  continuous  so  a 
small  error  in  the  coefficients  may  cause  a  large  fluctuation  in  the  representation.  Continuous  and 
n  times  differentiate  wavelets  are  known  [2]. 

The  two  equivalences  in  (3.8)  are  valid  for  finite-energy  deterministic  functions  in 
L2(R);  they  also  make  sense  for  finite-power  deterministic  function  in  L2(I),  compact  ICE, 
and  for  finite-energy  and  finite-power  random  processes  w.p.l  on  compact  ICR  when  compact 
scaling  function  4>(t)  and  compact  wavelet  function  y(t)  are  employed  [7],  [12]. 

3.3      Wavelet  Conditions  and  Generations 

It  is  clear  that  one  needs  to  obtain  the  wavelet  filter  sequences  h  and  g  in  order  to  perform 
the  D WT.    For  the  computation  of  the  CWT  in  (2. 1 ),  however,  one  needs  to  know  the  wavelet 
function  vj/(t)  explicitly.  For  compactly  supported  wavelets,  knowing  h  and  g  can  produce  the 
time-functions  <}>(t)  and  \j/(t).  In  this  section,  we  discuss  the  conditions  for  generating  an 
orthonormal  scaling/wavelet  pair,  and  give  details  for  how  to  produce  the  time  functions.  This 
section  is  important  for  the  fact  that  it  provides  working  knowledge  for  both  the  DWT  and  the 
CWT. 

3.3.1      Wavelet  conditions 

Consider  the  nested  sequence  of  closed  subspaces  { Vm;  m  e  2 }  as  discussed  in  section 
3.1.    The  scaling  function  (at  level  m)  <J>m(t)  is  in  Vm,  which  in  turn  is  a  proper  subset  of  the 
closed  space  Vm+i.  The  basis  of  Vm+i  is  given  by  {<J>m+i,n(t);  n  e  2 }  =  {2(m+l)/2<|>(2m+1 1  -  n); 
n  e  Z },  where  we  have  set  T  =  1  WLOG.  Since  4>m(t)  is  also  in  Vm+i,  we  can  represent  it  in 
terms  of  the  basis  functions  of  Vm+i: 

4>m(t)  =  X    <<t>m(0,  4>m+l,n(t)>    4>m+l,n(t)  (3-10) 


where  the  inner  product  <  <t>m(t)>  <t>m+l,n(t)  >  is  given  by 
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$m(t),<|>m+ifn(t)>=  J<Mt)  4>*m+l,n(t)   dt 


=  J  2m/2(j)(2mt  )  2(m+l)/2(})*(2m+lt  _n)  <jt 

=  V2  J  <f)(t)  4>(2t-n)  dt  (3.11) 

for  real  basic  scaling  function  <})(t).  Equation  (3.1 1)  describes  exactly  how  two  consecutive 
subspaces  Vm  C  Vm+i  are  correlated.  Note  that  (3.1 1)  depends  only  on  the  scaling  function  4>(t), 
but  not  on  the  individual  levels  m  and  m+1;  it  holds  for  any  value  of  m.  By  defining  h(n)  via  the 
following  dilation  equation: 

h(n)  =  J  <|>(t)  4>(2t  -n)   dt  (3.12) 

for  all  n,  we  can  rewrite  (3.10)  as 

4>m(t)  =  V2  £h(n)<|>m+i,n(t)  (3.13) 

n 

The  dilation  equation  of  (3. 12)  is  the  most  fundamental  equation  in  wavelet  theory.  One  can 
proceed  along  this  line  of  argument  and  derive  the  D WT  of  Figure  3.3.  For  a  scaling  function  <J>(t) 
that  has  a  finite  support,  the  number  of  solutions  to  (3.12)  is  finite,  and  so  h  is  of  finite  length. 

It  can  be  shown  that  [16]: 

J  <J)(t)  dt  =  constant  (3.14) 

and  WLOG  we  can  take  the  constant  to  be  1 .  This  implies  that  <D(f)  =  1  at  f  =  0,  which  means  that 
the  scaling  function  is  an  impulse  response  of  a  lowpass  type  filter.  Integrating  both  sides  of 
(3.13)  with  respect  to  (w.r.t.)  time  t  and  letting  m  =  0  yields: 


1  =  V2   Ih(n)    J2l/2<t>(2t-n)dt 
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=  X  h(n)  (3.15) 

n 

Equation  (3.15)  is  the  first  wavelet  condition,  which  is  a  consequence  of  the  nested  nature  of  the 

increasing  subspace  sequence  { Vm}.    For  orthonormal  scaling  function  <))(t),  (<})(t-nT);  n  e  Z }  is 

an  orthonormal  basis  for  Vo,  and  therefore, 

8(n)  =  <<Kt),<t>(t-n)> 

=  <  2  X  h(k)  4>(2t-k),   2  5>G-2n)  4>(2t-j)   > 
k  J 


-412 


h(k)hG-2n)55(k-j) 
k      j  2 


=2S 


h(k)h(k-2n)  (3.16) 

k 


Equation  (3.16)  is  the  second  wavelet  condition;  it  must  hold  if  the  scaling/wavelet  system  is 
orthonormal.  Performing  a  discrete  Fourier  transform  on  both  sides  of  (3. 1 6)  yields  the  conjugate 
filter  equation  of  (3.5).  Note  that  as  a  corollary  to  (3.16),  we  get  by  setting  n  =  0  that 


1 


=  X    h2(k)  (3.17) 


One  needs  only  to  find  h  such  that  it  satisfies  the  two  wavelet  conditions  of  (3.15)  and  (3.16). 
Once  h  is  known,  the  sequence  g  is  given  by  g(n)  =  (-l)1_n  h(l-n).  These  two  filter  sequences 
are  all  that  are  required  to  perform  the  DWT. 

3.3.2  Generation  of  wavelet  function 

In  the  previous  section,  we  have  shown  that  if  the  scaling  function  <j)(t)  is  known,  then  we 
can  determine  h  and  g  and  perform  the  DWT  with  ease.  Conversely,  if  h  is  found  by  solving  the 
two  wavelet  conditions,  then  we  can  generate  <J>(t)  and  y(t)  recursively.  The  construction  is  as 
follows.  Using  (3.13)  with  m  =  0,  we  get  the  fundamental  recursive  equation  for  the  scaling 
function  <}>(t): 


=  2l 


<t>(t)  =  2^  h(n)  <J)(2t-n)  (3.18) 

n 


Assuming  that  h  =  [h(0),  h(l),  h(2),  h(3), ,  h(M-l)],  i.e.,  the  h  sequence  has  finite  length, 

so  that  0(t)  has  a  compact  support  on  [0,  M-l],  i.e.,  0(t)  is  non-zero  on  this  compact  interval.  We 
first  find  the  values  of  0(t)  on  the  integers;  this  is  done  by  using  (3.18)  to  set  up  the  following 
matrix  equation: 


0(0) 
0(1) 

0(2) 


1_0(M-1)  J 


=  2 


h(0) 

0 

0 

0 

0 

h(2) 

h(l) 

h(0) 

0 

0 

h(4) 

h(3) 

h(2) 

h(l) 

h(0) 

• 

h(5) 

h(4) 

h(3) 

h(2) 

• 

• 

• 

h(5) 

h(4) 

h(M-2)  .... 
0       h(M-l)h(M-2)    • 
0  0  0       h(M-l)h(M-2) 

0  0  0  0  0 

0  0  0  0  0 

L0  0  0  0  0 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Ml) 

h(0) 

0 

h(3) 

h(2) 

h(l) 

h(5) 

h(4) 

h(3) 

h(5) 


h(M-l)h(M-2)    • 
0         0      h(M-l) 


0(0) 

0(D 

0(2) 


L<J>(M-1). 


for  even  M,  and 


(3.19a) 
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<j>(0) 

0(1) 

<l>(2) 


L<j)(M-l)  J 


=  2 


h(0) 

0 

0 

0 

0 

0 

h(2) 

h(l) 

h(0) 

0 

0 

0 

h(4) 

h(3) 

h(2) 

h(l) 

h(0) 

0 

• 

h(5) 

h(4) 

h(3) 

h(2) 

h(l) 

• 

• 

• 

h(5) 

h(4) 

h(3) 

h(M-l)h(M-2)    • 

0  0    h(M-l)  h(M-2)    • 

0  0  0  0    h(M-l)  h(M-2) 

0  0  0  0  0* 

0  0  0  0  0- 

L0  0  0  0  0 


0  0  0 

0  0  0 

0  0  0 

0  0  0 

h(0)  0  0 

h(2)  h(l)  h(0) 

h(4)     h(3)  h(2) 
h(4) 


h(M-l)h(M-2)    • 
0         0      h(M-l) 


<K0) 
0(2) 


L4>(M-1)  J 


(3.19b) 


for  odd  M.     In  general,  solving  (3.19a)  or  (3.19b)  is  equivalent  to  finding  the  eigenvector 
corresponding  to  the  eigenvalue  X  =  1  of  the  matrix  A  given  in  below: 

$  =  A$ 


(3.20) 


We  may  impose  the  additional  condition  on  <t>  that  ]T  <J)(n)  =  1 .  When  the  matrix  A  has  zero 

n 

determinant  there  can  be  multiple  X  =  1  solutions  to  (3.20).  Once  we  know  the  values  of  <j)(t)  on 
the  integers,  we  can  compute  its  values  on  the  ^-integers  by  using  the  fundamental  recursion  of 

(3.18)  as  follows: 


<K 


h  =  2Z,  h(n 


)  4>(k-n) 


(3.21) 


and  so-forth  for  on  the  t,  g,  -77,  ^  ...  -integers,  recursively.  Thus,  the  scaling  function  4>(t), 

which  has  support  on  the  compact  time  interval  [0,  M-l],  can  be  determined  explicitly  for  t  being  a 
dyadic  number.  Since  the  dyadic  numbers  (i.e.,  multiples  of  2~k,  for  all  k  >  0)  are  dense  in  R,  the 
function  <{>(t)  can  be  determined  at  all  times  t  in  the  limit. 
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The  wavelet  function  \|/(t)  is  similarly  computed  as  follows.    Since  \j/(t)  e  Wo  C  Vj,  we 
can  represent  \|/(t)  as  a  linear  combination  of  the  basis  functions  of  Vi : 

V(0  =  X  <  V(0.  21/2  <t>(2t-n)  >  2]/2  (J,(2t-n) 
n 

=  2X  g(n)4>(2t-n)  (3.22) 


which  is  a  recursion  similar  to  (3.18).  Since  we  already  know  <J)(t)  on  the  dyadic  numbers,  it 
follows  that  \|/(t)  can  be  computed  using  (3.22).  The  filter  sequence  g  is  a  flip-alternate-negation 

of  h,  and  we  can  set  g  =  [g(0),  g(l),  g(2), ,  g(M-l)]  so  that  \j/(t)  also  has  a  compact 

support  on  the  time  interval  [0,  M-l]. 

At  this  point  we  have  shown  how  to  generate  \j/(t)  for  the  CWT.  The  question  is:  How  do 
we  find  h(n)  in  the  first  place?   Although  one  may  guess  or  even  trial  and  error  to  find  a  h 
sequence  that  satisfy  the  two  wavelet  conditions,  in  reality  it  is  a  difficult  task  to  find  a 
scaling/wavelet  pair  that  has  nice  analytical  properties.   Nonetheless,  several  researchers  have 
succeeded  in  systematically  generating  wavelet  functions.  We  show  some  of  their  work  in  the 
following. 

3.3.3  Compactly  supported  wavelets  examples 

First  of  all,  the  simplest  compactly  supported  wavelet  is  the  Haar  wavelet.  It  has  a  compact 
support  on  [0,  1],  with  h  =  [1/2,  1/2]  and  g  =  [1/2,  -1/2].  The  Haar  scaling/wavelet  pair  is  shown 
in  Figure  3. 4. a.  Note  that  the  functions  are  not  continuous. 

Daubechies  [2]  has  proposed  a  class  of  compactly  supported  orthonormal  wavelets  and  she 
gave  the  hs  for  M  =  4,6,8,10,12,14,16,18,20.  For  example,  the  Daubechies  4-coefficient  wavelet 
has  h  =  [h(0),  h(l),  h(2),  h(3)]  given  by: 

uffvs  _  1+^3  .  3+V3 

h(0)  =  — g—   ,        n(l)=  — g— 

3-V3  1-V3 

h(2)  =  ^jp   ,         h(3)  =  -^  (3.23) 
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Note  that  h(0)+h(l)+h(2)+h(3)  =  1,  and  h2(0)+h2(l)+h2(2)+h2(3)  =  1/2.  Substituting  h  into 
(3.19)  we  find  the  eigenvector  for  X  -  1  to  be 

M>(0),  <K1),  4>(2),  4>(3)]  =  [0,  1.366,  -0.366,  0]  (3.24) 

The  corresponding  scaling  function  <)>(t)  and  the  wavelet  function  \y(t)  (D4)  are  recursively 
computed,  using  (3.18)  and  (3.22),  respectively,  and  they  are  plotted  in  Figure  3.4.b  We  also 
show  the  Daubechies  10-coefficient  scaling/wavelet  (D10)  pair  in  Figure  3.4.c  and  the  20- 
coefficient  scaling/wavelet  (D20)  pair  in  Figure  3.4.d  Observe  that  D20  looks  like  a  modulated 
Gaussian.  Also  note  that  the  Daubechies  wavelets  are  increasingly  smoother  in  time  and  better- 
contained  in  the  frequency  spectrum. 

Pollen  [19]  has  shown  that  it  is  possible  to  parametrize  compactly  supported  wavelets.  In 
particular,  for  M  =  6,  using  a  pinched  torus,  he  was  able  to  give  a  parametriztion  of  the  h  sequence 
as  follows: 

h(0)  =  [(1  +  cos  a  +  sin  a)(l  -  cos  (3  -  sin  p)  +  2  sin  p  cos  a  ]/8 

h(l)  =  [(1  -  cos  a  +  sin  a)(l  +  cos  p  -  sin  p)  -  2  sin  p  cos  a  ]/8 

h(2)  =  [1  +  cos  (a-p)  +  sin  (a-p)]/4 

h(3)  =  [1  +  cos  (a-P)  -  sin  (a-p)]/4 

h(4)  =  [l-h(0)-h(2)]/2 

h(5)  =  [l-h(l)-h(3)]/2  (3.25) 

where  -n  <  a,  P  <  n  are  the  two  parametrizing  angles.  It  is  indeed  a  very  nice  and  important 
result.    As  an  example,  the  Haar  wavelet  has  parametrization  a  =  P,  for  any  a.  As  another 
example,  we  randomly  generate  a  =  0.2172  and  P  =  0.3359  and  use  (3.25)  to  obtain 

h  =  [0.0055,  -0.0322,  0.4686,  0.5278,  0.0259,  0.0044].  (3.26) 

It  turns  out  that  this  is  a  valid  wavelet  sequence;  we  carry  out  the  recursive  generation  and  we  plot 
the  resulting  scaling  function  and  wavelet  function  in  Figure  3.4.e.  This  wavelet  has  dual 
passbands.  Note  that  the  scaling  function  and  the  wavelet  function  are  orthonormal  basis 
functions.  One  can  experiment  with  (3.26)  to  generate  very  crazy  looking  wavelet  functions. 
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Figure  3.4.a     The  Haar  scaling  function  and  wavelet  function. 
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Figure  3.4.b     Daubechies  4-coefficient  (D4)  scaling  function  and  wavelet  function, 
and  their  energy  spectra. 
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Figure  3.4.C     Daubechies  10-coefficient  (Din)  scaling  function  and  wavelet  function, 
and  their  energy  spectra. 
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Figure  3.4.d    Daubechies  20-coefficient  (D20)  scaling  function  and  wavelet  function, 
and  their  energy  spectra. 


It 


phi(t) 


PSD  phi(t) 


PSD  xi(t) 


time 


-2  0  2 

frequency 


Figure  3.4.e     A  6-coefficient  scaling  function  and  wavelet  function  with  Pollen 
parametrization  a  =  0.2172,  P  =  0.3359,  and  their  energy  spectra. 
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3.4      Multiple-Phase  Discrete  Wavelet  Transform 

In  digital  signal  analysis,  what  one  has  at  hand  is  a  sequence  of  N  samples,  x(nT),  where  T 
is  the  sampling  period  and  n  =  0,  1,2,...  Ns-1.  It  is  convenient  to  consider  an  auxiliary  function 
xo(t)  in  Vo:  xo(t)  =  J^  bu.n  to.nW  by  setting  bo,n  =  x(nT)  at  resolution  level  0,  which  is  the 

n 

highest  level  available.  Using  a  recursive  filtering  algorithm,  the  so-called  discrete  wavelet 
transform  (DWT)  of  section  3.2,  one  can  compute  the  lower-level  resolution  coefficients  as  (Fig. 
3.3): 

bm-i(i)=  X  >/  2  h(k)  bm(2i  +  k) 

k 

dm.i(i)  =  X  V  2  g(k)  bm(2i  +  k)  (3.27) 

k 

for  m  =  0,  -1,  -2,  -3 Here,  h  =  [  h(-M+2),  h(-M+3) ,  h(0),  h(l)  ]  and  g  =  [  g(-M+2), 

g(-M+3), . . . ,  g(0),  g(l)  ]  are  the  appropriate  conjugate-quadrature  discrete-time  filter  responses. 

oo 

The  h  coefficients  are  obtained  by  the  dilation  equation  h(k)  =    J  <}>(t)  <J)(2t  -  k)  dt.    We  have  set 

-oo 

g(k)  =  (-i)-M+3-k  h(-M+3-k)  so  that  the  time  indexes  for  the  b's  and  d's  are  always  non-negative. 
We  also  assumed  that  both  $  and  y  have  compact  supports  so  that  h  and  g  have  only  finite  number 
of  filter  taps  (  =  M).  Roughly  speaking,  the  approximation  coefficients  bm,n  capture  the 
approximate  information  in  the  signal  sequence  at  resolution  level  m,  while  the  difference 
coefficients  dm>n  retain  the  detail  information  going  from  resolution  level  m-  \  torn.  The  DWT 
has  important  applications  in  signal  analysis,  particularly  in  compression,  detection,  and 
classification.  The  choice  of  <{)  and  \y  (hence  h  and  g)  is  such  that  the  signal  is  localized  as  much  as 
possible  in  both  time  and  frequency.  Unfortunately,  the  DWT  is  not  time  invariant.  If  the  signal 

sequence  is  time-shifted  by  one  unit  to  become  t>  =  [  0,  bo,o»  bo.i,  bo,2» bo,N-l  ] tnen  tne 

corresponding  DWT  is  not  time-shifted  and  it  will  be  completely  different  from  the  old  DWT.  This 
will  not  present  a  problem  to  the  reconstruction  of  the  signal  sequence;  however,  the  merits  of 
DWT  for  signal  detection  and  classification  may  be  diminished. 

We  note  that  for  dyadic  DWT,  there  are  two  different  sets  of  DWT  at  each  level;  we  name 
them  as  "phase-0"  and  "phase- 1"  decompositions.  Referring  to  Fig.  3.5,  for  phase-0 
decomposition,  the  0-th  coefficient  at  resolution  level  m  -  1  is  [h(l)]  •  [bm)U],  me  l~st  coefficient  is 
[h(-l),  h(0),  h(l)]  •  [bm,o,  bm,i,  bm,2l  where  "•"  denotes  dot-product,  the  2-nd  coefficient  is 
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[h(2),  h(-l),  h(0),  h(l)]  •  [bm,i,  bm,2,  bm,3,  bm,4],  and  so  forth.  On  the  other  hand,  the  phase- 1 
coefficients  will  read  as  [h(0),  h(l)]  •  [bm,n,  bm,i],  [h(-2),  h(-l),  h(0),  h(l)]  •  [bm,0,  bm,i,  bm,2, 
bm,3L  [h(-2)t  h(-l),  h(0),  h(l)]  •  [bm,2,  0^3, 1^,4,  0^5] 
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Figure  3.5       Two  phases  of  wavelet  decomposition. 


We  have  proposed  a  multiple-phase  DWT  (MP/DWT)  for  signal  analysis  [10],  [14],  [13].  The  key 
to  this  MP/DWT  is  to  construct  new  filter  sequences  hm  and  gm  (m  <  -1)  at  each  decomposition 
level  by  inserting  2*(m+I)  -  1  number  of  zeros  between  every  two  subsequent  filter  coefficients. 
The  MP  approximation  coefficients  at  level  m  -  1,  time  indexed  i,  is  given  by  the  dot  product  of  hm 

1  1  1 

and  b  ,  where  b  is  the  MP  approximation  coefficient  sequence  at  level  m  denoted  by  b  = 

1         1         1 
(bm.o-  bm,r  bm.2  • >'  Specifically,  we  have 


m-l,i 


1 

I 


V  2  hm(k)  b 


m,i+k-l 


and 


k=  -M+2-(M-l)(2-m-l) 


1 


Cu  =     x v  2  gm(k)  b"^k-1 

k=  -M+2-(M-l)(2-m-l) 


(3.28) 


The  manner  that  the  phase-0  and  phase- 1  coefficients  are  interleaved  at  each  level  may  be  described 
by  the  tree  diagram  shown  in  Fig.  3.6.   For  a  k-level  decomposition,  there  are  altogether  2k 
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possible  phases  of  decomposition.  Each  path  in  the  tree  of  Fig.  3.6  represents  an  unique  phase- 
sequence.  One  may  choose  to  decompose  the  signal  along  any  one  of  these  2k  phase-sequences 
and  obtain  a  different  DWT.  Interleaving  all  of  them  yields  our  MP/DWT.  Note  that  knowing  the 
phase-sequence,  one  can  trace  back  in  the  tree  and  reconstruct  the  signal  from  the  DWT.   Note  that 
MP/DWT  is  a  redundant  transformation.  It  may  be  regarded  as  a  discrete  approximation  of  the 
continuous  wavelet  transform  for  xn(t).  Its  time-invariant  property  may  be  useful  for  some  signal 
detection  problems.  The  MP/DWT  can  be  used  to  detect  PSK  signals,  but  we  shall  report  the 
results  elsewhere.  For  data  compression  applications,  one  can  maximize  the  data  compression 
ratio  by  optimizing  over  the  2k  phase-sequences. 
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Figure  3.6     Different  phases  of  wavelet  decomposition. 
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4 .        DETECTION  OF  PSK  SIGNALS 

We  now  use  the  CWT  of  (2.1)  to  analyze  PSK  signals.  Viewing  the  CWT  in  a  time-scale 
plane,  we  attempt  to  identify  qualitative  and  quantitative  features  that  would  reveal  the  exact 
occasions  of  the  phase  shifts.  First,  we  consider  a  basebad  rectangular  wave  and  its  noisy 
versions.  Next,  we  analyze  a  180°  phase  shift  in  a  sinusoid  and  its  noisy  versions.  Haar, 
Daubechies,  and  complex  modulated  Gaussian  wavelets  will  be  employed  in  the  decompositions. 

4.1      Rectangular  Wave 

4.1.1  Haar  wavelet 

A  baseband  rectangular  wave  and  its  CWTs  using  the  Haar  wavelet  are  shown  in  Figure 
4.1. a.   The  rectangular  wave  signal  has  a  normalized  duration  of  1  second.  The  signal  sequence 
x(nT)  has  length  equal  to  65  and  it  changes  sign  at  n  =  9,  25, 41,  and  57,  which  corresponds  to  t  = 
nT  =  9/64,  25/64, 41/64  and  57/64  seconds,  respectively.  The  Haar  wavelet  \j/(t)  =  \j/o(t)  was 
shown  in  Figure  4.  La.  The  transforms  CWx(nT,  2_m)  is  computed  using  (2.1)  for  n  =  0,1, .  . .  . 
.  ,64,  and  m  =  0,  1,  2, .  ...  ,6;  they  are  shown  in  sequential  order  in  Figure  4.1. a.  Note  that 
there  are  some  end  effects  due  to  the  finite  signal  window  and  we  would  ignore  these  end  effects. 
At  low  level  of  decomposition,  i.e.,  m  =  0,1,  or  2,  the  CWT  has  good  frequency  resolution  and 
we  see  that  the  oscillating  frequency  of  the  signal  is  revealed.    At  higher  level  of  decomposition, 
the  time  resolution  is  getting  better  and  the  CWT  reveals  the  sign  change  instances  with  increasing 
accuracy.    Specifically,  at  m  =  6,  a  positive  "spike"  in  the  CWT  shows  a  +1  to  -1  sign  change, 
and  a  negative  "spike"  shows  a  -1  to  +1  sign  change.  The  Haar  wavelet  is  matched  to  the  sign 
change  so  the  good  time  resolution  is  expected.  However,  at  high  time  resolution,  the  wavelet 
filter  has  large  bandwidth.   When  there  is  noise,  a  lot  of  noise  energy  will  be  collected  and  the 
signal  CWT  will  be  smeared,  as  we  will  show  in  the  next  figure. 

We  add  noise  to  the  rectangular  wave  and  show  the  CWTs  in  Figure  4. 1  .b  and  4. 1  .c. 
Power  of  the  signal  sequence  is  E(x2(nT))  =  1.  The  random  noise  sequence  is  Gaussian  with 
mean  zero  and  variance  s2.  Define  the  signal-power  to  noise-power  ratio  as: 

SNR  =    signal  sequence  power 

noise  sequence  power  v  '  ' 
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In  Figure  4.1.b,  the  SNR  is  10  dB  and  in  Figure  4.1.C  the  SNR  is  5  dB.  We  see  that  the  low  level 
CWTs  are  not  severely  affected  so  we  still  have  good  frequency  resolution.  At  higher  level,  the 
CWTs  are  increasingly  smeared  by  the  noise,  especially  at  m  =  6.  This  is  because  the  output 
signal-power  to  noise-power  ratio,  SNRq,  defined  by: 

<;isfR     _  filtered  signal  sequence  power  ..  «. 

^  "~    filtered  noise  sequence  power 

is  decreasing  exponentially  with  m.   Note  that  the  wavelet  filter  (with  impulse  response  Vm(0)  has 
bandwidth  B  =  2»2m  Go  and  so  the  filtered  noise  power  is  increasing  exponentially  with  m  for 
wideband  noise.  The  filtered  signal  power,  however,  will  decrease  with  m  since  the  r.m.s.  center 
frequency  of  the  wavelet  filter  is  moving  towards  f  =  -H»  and  smaller  and  smaller  amount  of  signal 
energy  will  be  collected.  If  the  signal  is  a  wideband  spread-spectrum  signal  instead  of  a 
narrowband  signal,  then  it  is  possible  that  SNRo  will  remain  large  for  many  levels  of 
decompositions.  This  suggests  that  WT  is  suitable  for  detecting  spread-spectrum  signals. 

For  SNR  =  10  dB  as  in  Figure  4. 1  .b,  one  can  observe  from  the  CWT  at  m  =  4  and  m  =6 
that  there  are  4  sign  changes.  However,  when  one  attempts  to  pinpoint  the  exact  times  of 
occurrence  by  viewing  the  CWT  at  m  =  6,  only  the  changes  at  t  =  9/64  and  41/64  can  be  located 
with  confidence. 

For  SNR  =  5  dB  in  Figure  4. 1  .c,  the  noise  power  is  too  much  for  the  CWT  to  reveal  good 
time-resolutions.  Only  the  oscillating  frequency  of  the  rectangular  wave  can  be  revealed.  The 
Haar  wavelet  does  not  work  well  in  noise  because  its  poor  spectral  properties,  i.e.,  its  spectrum 
decays  too  slow  and  has  too  many  sidelobes. 
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Figure  4. 1  .a    CWT  of  rectangular  wave  with  Haar  wavelet.    No  noise. 
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Figure  4.1.b    CWT  of  rectangular  wave  with  Haar  wavelet.    SNR  =  10  dB. 
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Figure  4. 1  .c    CWT  of  rectangular  wave  with  Haar  wavelet.    SNR  =  5  dB. 
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4.1.2  Daubechies  wavelet 

Daubechies  wavelets  have  better  frequency  properties  than  the  Haar  wavelet.    However, 
they  are  not  matched  (in  time)  to  the  sign  change  in  the  rectangular  wave.  We  show  the  CWTs  for 
the  rectangular  wave  with  Daubechies  10-coefficient  wavelet  (Dio)  in  Figures  4.2.a,  4.2.b,  and 
4.2.C  for  the  noiseless,  SNR  =  10  dB,  and  SNR  =  5  dB  cases,  respectively.  The  wavelet  is  first 
normalized  in  time  so  that  it  is  confined  to  the  interval  [0, 1]  at  level  m  =  0. 

When  there  is  no  noise,  the  Daubechies  wavelet  is  also  capable  in  locating  the  sign  changes 
with  high  time-resolution  as  shown  in  Figure  4.2. a.  Furthermore,  because  its  frequency  spectrum 
is  more  confined  than  the  Haar  wavelet,  the  CWT  at  low  level  reveals  the  oscillating  frequency  of 
the  rectangular  wave  with  greater  resolution. 

When  SNR  =  10  dB,  the  performance  is  worse,  as  shown  in  Figure  4.2.b.  At  levels  m  =  3 
and  m  =  4,  it  is  possible  to  identify  the  first  and  the  third  sign  changes  correctly,  just  as  in  the  Haar 
case.  But  there  is  a  false  indication  around  t  =  0.4.    When  SNR  =  5dB,  this  Daubechies  wavelet 
is  seen  to  perform  not  much  better  than  the  Haar  wavelet,  as  shown  in  Figure  4.2.c. 

The  problem  with  both  Haar  and  Daubechies,  as  well  as  other  unmodulated  wavelets  is  that 
they  are  not  turnable.  It  is  difficult  to  select  or  "customize"  the  passband  center  frequency  for  a 
specific  signal,  thereby  the  wavelet  filter  may  miss  most  of  the  signal  energy.  It  is  also  not 
possible  to  have  a  relatively  narrow  bandwidth  at  high  frequencies.  The  complex  modulated 
Gaussian  wavelet  is  capable  to  alleviate  these  shortcomings. 
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Figure  4.2.a    CWT  of  rectangular  wave  with  Daubechies  10-coefficient  wavelet. 
No  noise. 
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Figure  4.2.b    CWT  of  rectangular  wave  with  Daubechies  10-coefficient  wavelet. 
SNR=  10  dB. 
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Figure  4.2.c    CWT  of  rectangular  wave  with  Daubechies  10-coefficient  wavelet. 
SNR  =  5dB. 
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4.1.3  Complex  modulated  Gaussian  wavelet 

By  changing  the  modulating  frequency,  fy,  of  the  complex  modulated  Gaussian  wavelet, 
one  can  place  the  wavelet  passband  center  frequency  practically  anywhere  on  the  positive 
frequency  axis.  Changing  fy,  as  well  as  a,  in  (2.10)  gives  us  two  extra  dimensions  to  expose  the 
signal  features.  Although  there  may  be  a  lot  of  redundancy,  this  approach  essentially  gives  us  a 
tunable,  variable-bandwidth  multiresolutional  analysis  of  the  signal.  It  is  also  possible  to  adapt  the 
analysis  to  different  types  of  signal  and  noise. 

For  simpler  presentation,  we  fix  a  =  0.12  and  only  vary  fy  to  be  4,  8,  and  16  Hz.  For  the 
noiseless  case,  we  show  the  corresponding  time-scale  decompositions  (which  are  complex  valued) 
in  Figures  4. 3. a  -  4.3. f.  The  first  three  figures  depict  the  magnitudes  of  the  CWTs,  while  the  last 
three  figures  show  the  phases  of  the  CWTs.    Consider  the  magnitude  plots  in  Figures  4. 3. a,  4.3.b, 
and  4.3.C.  There  are  altogether  3x7  =  21  snapshots  of  the  wavelet  filtered  rectangular  wave.  Each 
set  of  7  snapshots  start  with  a  different  center  frequency  fv  and  the  filter  banks  scan  the  frequency 
axis  as  m  increases  from  0  to  6.  Most  of  these  21  magnitude  plots  do  correctly  reveal  the  sign 
change  times  (with  different  time-resolutions)  and  would  be  useful  in  the  noisy  cases.  The  phase 
plots  in  Figures  4.3.d,  4.3.e,  and  4. 3. f  confirm  the  sign  change  times,  and  in  addition  they  also 
reveal  the  polarity  (+  to  -  or  -  to  +)  of  each  change.  Note  that  when  there  is  noise  and  the 
magnitude  plots  are  not  conclusive,  we  can  consult  the  phase  plots  for  a  confirmation.   This 
alternative  source  of  information  is  not  available  in  decompositions  using  real  wavelets. 

For  the  case  when  SNR  =  10  dB,  we  show  the  magnitude  plots  for  fy  =  4,  8,  16  in 
Figures  4.4.a,  4.4.b,  and  4.4.c.  In  Figure  4.4. b  at  m  =  6,  the  magnitude  plot  of  the  CWT  clearly 
indicate  that  there  are  sign  changes  at  the  four  correct  instance,  although  it  addition  there  is  a 
possible  false  indication  at  t  *  0.4.  Figure  4.4.a  at  m  =  3, 4,  and  5  confirm  that  the  findings  are 
due  to  the  signal  but  not  the  noise.  Note  that  in  Figure  4.4.c,  the  magnitude  plots  are  not 
informative  and  certainly  are  not  conclusive;  the  choice  of  fy  =  16  is  too  large  for  this  particular 
signal. 

When  SNR  =  5  dB,  we  show  one  set  of  magnitude  plots  of  the  CWT  in  Figure  4.5  with  fy 
=  8.    At  level  m  =  6,  the  CWT  reveals  correctly  the  first,  second,  and  the  third  sign  change, 
probably  misses  the  forth  one  and  falsely  add  one  at  t  =  0.4.  The  results  can  then  be  referenced 
with  other  magnitude  plots  using  a  different  modulating  frequency.  Comparing  with  the  Haar  and 
the  Daubechies  wavelets,  the  performance  of  the  complex  modulated  Gaussian  is  better,  especially 
in  the  noisy  cases. 
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Figure  4.3.a    Magnitude  of  CWT  of  rectangular  wave  with  complex  modulated 
Gaussian  wavelet.  o  =  0.12,  fy  =  4.  No  noise. 
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Figure  4.3. b    Magnitude  of  CWT  of  rectangular  wave  with  complex  modulated 
Gaussian  wavelet.  o  =  0.12,  fy  =  8.  No  noise. 
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Figure  4.3.C    Magnitude  of  CWT  of  rectangular  wave  with  complex  modulated 
Gaussian  wavelet,  a  =  0.12,  fv  =  16.  No  noise. 
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Figure  4.3.d    Phase  of  CWT  of  rectangular  wave  with  complex  modulated 
Gaussian  wavelet,  a  =  0.12,  fy  =  4.  No  noise. 
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Figure  4.3.e    Phase  of  CWT  of  rectangular  wave  with  complex  modulated 
Gaussian  wavelet,  a  =  0.12,  fv  =  8.  No  noise. 


45 


signal 


CWT,  level  2A-3 


200 


-200 


200 


-200 


200 


0.5 

time 

CWT,  level  2^ 


0.5 
CWT,  level  2A-1 


0.5 
CWT,  level  2A-2 


200 


-200 


200 


-200 


200 


-200 


200 


-200 


0.5 
CWT,  level  2A-4 


0.5 
CWT,  level  2A-5 


0.5 
CWT,  level  2A-6 


r 


0.5 


Figure  4.3.f    Phase  of  CWT  of  rectangular  wave  with  complex  modulated 
Gaussian  wavelet,  a  =  0.12,  fv=  16.  No  noise. 
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Figure  4.4.a    Magnitude  of  CWT  of  rectangular  wave  with  complex  modulated 
Gaussian  wavelet,  a  =  0.12,  fy  =  4.    SNR  =  10  dB. 
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Figure  4.4.b    Magnitude  of  CWT  of  rectangular  wave  with  complex  modulated 
Gaussian  wavelet,  a  =  0.12,  fv  =  8.    SNR=10dB. 
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Figure  4.4.c    Magnitude  of  CWT  of  rectangular  wave  with  complex  modulated 
Gaussian  wavelet,  a  =  0.12,  fv  =  16.    SNR  =  10  dB. 
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Figure  4.5       Magnitude  of  CWT  of  rectangular  wave  with  complex  modulated 
Gaussian  wavelet.  a  =  0.12,  fv  =  8.    SNR  =  5  dB. 
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4.2      180°     phase  shift 

The  signal  to  be  analyzed  is  now  changed  to  a  sinusoid  with  a  1 80°  phase  shift  at  t  =  0,  as 
shown  in  the  first  subplot  of  Figure  4.6.  We  use  the  Haar,  Daubechies  10-coefficient,  and  the 
complex  modulated  Gaussian  wavelets  to  analyze  this  signal  and  plot  their  CWTs  in  Figures  4.6, 
4.7,  and  4.8,  respectively.  For  the  Haar  wavelet,  the  phase  change  also  appears  as  a  phase  change 
in  most  of  the  CWTs  as  shown  in  Figure  4.6.  The  Daubechies  wavelet  does  a  better  job  in 
identifying  the  phase  change,  as  shown  in  Figure  4.7.  The  complex  modulated  Gaussian  wavelet 
also  reveals  clearly  the  phase  shift  in  the  magnitude  plots  of  the  CWT  at  m  =  4  and  m  =5,  in 
addition  a  phase  change  is  confirmed  by  m  =  1  and  m  =2  plots,  as  shown  in  Figure  4.8. 

When  there  is  noise.,  we  found  that  both  the  Haar  and  Daubechies  wavelets  failed  to 
perform. 

When  SNR  =  10  dB  (note  that  the  signal  power  is  equal  to  0.5),  we  choose  a  =  0.12  and 
fy  =  12  for  the  complex  modulated  Gaussian  and  show  the  magnitude  and  phase  plots  in  Figures 
4.9. a  and  4.9.b.  In  Figure  4.9. a,  the  magnitude  plot  of  the  CWT  at  m  =  2  reveals  a  phase  shift 
activity  somewhere  on  or  before  t  =  0.5.  (Note  that  there  is  an  end-effect  in  the  DWT  due  to  finite 
signal  window.)   But  the  rest  of  the  magnitude  plots  fail  to  pinpoint  a  better  time  location.  When 
we  consult  the  phase  plots  as  shown  in  Figure  4.9.b,  the  m  =  6  plot  suggests  that  there  is  an 
abrupt  change  of  phase  at  t  =  0.5.  This  is  confirmed  by  the  m  =  3  and  m  =  6  phase  plots  in  Figure 
4.9.c,  where  we  use  a  different  f y  =  16. 

When  SNR  =  5  dB,  we  choose  a  =  0.12,  fy  =  8.  The  magnitude  subplot  for  m  =  2  of 
Figure  4.10.a  reveals  some  phase  activities;  however,  the  other  magnitude  subplots  and  the  phase 
plots  of  Figure  4.10.b  yield  no  confirmative  information.  Detection  of  a  single  phase  shift  is  very 
sensitive  to  the  noise.  However,  when  the  signal  window  is  much  longer  and  more  phase  shifts 
are  present,  the  complex  modulated  Gaussian  would  perform  due  to  a  larger  output  signal-power  to 
noise-power  ratio  (SNRo,  see  (4.2) ). 
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Figure  4.6     Magnitude  of  CWT  of  1 80°  phase  shifted  sinusoid  with  Haar 
wavelet.  No  noise. 
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Figure  4.7     Magnitude  of  CWT  of  1 80°  phase  shifted  sinusoid  with  Daubechies 
10-coefficient  wavelet.  No  noise. 
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Figure  4.8     Magnitude  of  CWT  of  1 80°  phase  shifted  sinusoid  with  complex 
modulated  Gaussian  wavelet,  o  =  0. 12,  fy  =  8.    No  noise. 
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Figure  4.9.a     Magnitude  of  CWT  of  180°  phase  shifted  sinusoid  with  complex 
modulated  Gaussian  wavelet,  a  =  0.12,  f¥  =  12.    SNR  =  10  dB. 
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Figure  4.9.b     Phase  of  CWT  of  1 80°  phase  shifted  sinusoid  with  complex 

modulated  Gaussian  wavelet,  c  =  0.12,  fv  =  12.    SNR  =  10  dB. 
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Figure  4.9.C     Phase  of  CWT  of  180°  phase  shifted  sinusoid  with  complex 

modulated  Gaussian  wavelet,  a  =  0.12,  fy  =  16.    SNR  =  10  dB. 
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Figure  4.  lO.a     Magnitude  of  CWT  of  1 80°  phase  shifted  sinusoid  with  complex 
modulated  Gaussian  wavelet,  a  =  0.12,  fv  =  8.    SNR  =  5dB. 
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Figure  4.  lO.b     Phase  of  CWT  of  1 80°  phase  shifted  sinusoid  with  complex 

modulated  Gaussian  wavelet,  a  =  0.12,  fy  =  8.    SNR  =  5  dB. 
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5.        CONCLUDING  REMARKS 

Wavelet  decompositions  have  better  time-resolution  at  higher  frequencies.  They  have 
better  frequency-resolution  at  lower  frequencies.  Wavelet  decompositions  are  suitable  for 
analyzing  signals  that  have  high-frequency  short  pulses  and/or  low-frequency  long  pulses.  The 
discrete  wavelet  transform  (DWT)  is  primarily  useful  for  data  compaction,  although  we  have  found 
some  of  its  potentials  in  signal  detection.  The  continuous  wavelet  transform  (CWT)  can  be  used  in 
signal  detection.  We  have  shown  that  by  using  complex  modulated  Gaussian  wavelets,  the  time 
and  frequency  features  of  a  signal  can  be  revealed  via  a  tunable,  variable-bandwidth, 
multiresolution  analysis.  We  have  shown  that  CWT  can  be  used  in  detecting  PSK  signals. 
However,  the  performance  is  sensitive  to  the  level  of  noise. 
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